Мягкий дедлайн: 4 ноября 2020 г. 23:59
Жесткий дедлайн (со штрафом в 50% от количества набранных вами за ДЗ баллов): 5 ноября 2020 г. 08:59
Визуализация "чего-либо" без выполненного основного задания оценивается в 0 баллов
ФИО: Гудков Дмитрий Всеволодович
Группа: DS-11
Генерируйте рандомные точки на планете Земля до тех пор, пока не попадете на территорию Афганистана
import nbconvert
jupyter nbconvert --to html HW2_.ipynb
File "<ipython-input-3-c97678e9646a>", line 1 jupyter nbconvert --to html HW2_.ipynb ^ SyntaxError: invalid syntax
from OSMPythonTools.nominatim import Nominatim
import overpy
import shapely.geometry as geometry
from shapely.ops import linemerge
import folium
import random
nominatim = Nominatim()
areaId = nominatim.query('Afganistan').areaId()
areaId
3600303427
query = """[out:json][timeout:25];
rel(303427);
out body;
>;
out skel qt; """
api = overpy.Overpass()
result = api.query(query)
lss = [] #convert ways to linstrings
for ii_w, way in enumerate(result.ways):
ls_coords = []
for node in way.nodes:
ls_coords.append((node.lon, node.lat)) # create a list of node coordinates
lss.append(geometry.LineString(ls_coords)) # create a LineString from coords
merged = linemerge([*lss]) # merge LineStrings
afganistan = geometry.Polygon(merged)
way = []
lat = random.uniform(-90, 90)
long = random.uniform(-180, 180)
way.append(geometry.Point(long, lat))
while True:
if afganistan.contains(geometry.Point(long, lat)):
break
if (lat <= afganistan.centroid.y) and (long <= afganistan.centroid.x):
lat = random.uniform(lat, afganistan.centroid.y)
long = random.uniform(long, afganistan.centroid.x)
elif (lat >= afganistan.centroid.y) and (long >= afganistan.centroid.x):
lat = random.uniform(afganistan.centroid.y, lat)
long = random.uniform(afganistan.centroid.x, long)
elif (lat <= afganistan.centroid.y) and (long >= afganistan.centroid.x):
lat = random.uniform(lat, afganistan.centroid.y)
long = random.uniform(afganistan.centroid.x, long)
elif (lat >= afganistan.centroid.y) and (long <= afganistan.centroid.x):
lat = random.uniform(afganistan.centroid.y, lat)
long = random.uniform(long, afganistan.centroid.x)
way.append(geometry.Point(long, lat))
map_of_afganistan = folium.Map([33.83133703361575, 66.02693526354392], zoom_start=2)
for point in range(len(way)):
folium.Marker(
location = [way[point].y, way[point].x],
popup = str(point) + ' step'
).add_to(map_of_afganistan)
folium.GeoJson(afganistan).add_to(map_of_afganistan)
map_of_afganistan
Визуализируйте пошагово предложенный алгоритм при помощи Folium
Для измерения показателя качества жизни в точке, найденной в предыдущем задании, вам необходимо рассчитать следующую сумму расстояний (в метрах):
* При поиске N ближайших объектов обязательно использовать R-tree
** Южной частью страны является территория, находящаяся к югу от множества точек, равноудаленных от самой северной и самой южной точек страны
from OSMPythonTools.overpass import Overpass
overpass = Overpass()
result = overpass.query('relation["admin_level"="2"][boundary=administrative];out;')
areas = []
for el in result.elements():
el_id = str(el.id())
if el_id[-1] == '0':
areas.append([el.id(), el.tag('name:en')])
from OSMPythonTools.overpass import overpassQueryBuilder, Overpass
counts = []
for area in areas:
areaId = 3600000000 + area[0]
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"landuse"="residential"', out='count')
result = overpass.query(query, timeout=900)
counts.append([result.countElements(), area[0], area[1]])
def maximum_residential(counts):
maximum = counts[0]
for country in counts:
if maximum[0] < country[0]:
maximum = country
return maximum
max_residential = maximum_residential(counts)
print(f'В {max_residential[2]} (id: {max_residential[1]}) содержится наибольшее количество объектов жилой недвижимости - {max_residential[0]}')
cameroon_id = max_residential[1]
В Cameroon (id: 192830) содержится наибольшее количество объектов жилой недвижимости - 139054
areaId = 3600000000 + cameroon_id
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"amenity"="atm"')
result = overpass.query(query, timeout=900)
points = []
for atm in result.elements():
points.append((atm.lon(), atm.lat()))
import rtree.index
idx = rtree.index.Rtree()
for i, p in enumerate(points):
idx.insert(i, p + p, p)
p = (way[-1].x, way[-1].y)
five_nearest = list(idx.nearest(p, 5, objects='raw'))
from shapely.geometry import Point, LineString
orig_point = Point(way[-1])
lines = [LineString([orig_point, Point(dest)]) for dest in points if dest in five_nearest]
from pyproj import Geod
geod = Geod(ellps="WGS84")
lines_length = [geod.geometry_length(line) for line in lines]
print(f'{lines_length} расстояния от точки до 5 ближайших банкоматов в Камеруне')
[6144389.378466479, 6557451.369484478, 6557127.224546507, 6514995.396643748, 6489205.27622738] расстояния от точки до 5 ближайших банкоматов в Камеруне
sum_length = 0
for i in lines_length:
sum_length += i
print(f'{sum_length} - сумма расстояний в метрах от точки до 5 ближайших банкоматов в Камеруне')
32263168.64536859 - сумма расстояний в метрах от точки до 5 ближайших банкоматов в Камеруне
countries_and_capitals = []
for area in areas:
areaId = 3600000000 + area[0]
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"capital"="yes"')
result = overpass.query(query, timeout=900)
if result.toJSON()['elements'] == []:
continue
else:
countries_and_capitals.append([area[0],
area[1],
result.toJSON()['elements'][0]['id'],
result.toJSON()['elements'][0]['tags']['name:en']])
countries_and_relation_capitals = []
for country in countries_and_capitals:
if country[1] == 'Brazil':
areaId = 3600000000 + 2758138
elif country[1] == 'Australia':
areaId = 3600000000 + 2354197
else:
nominatim = Nominatim()
capital_country = country[3] + ', ' + country[1]
areaId = nominatim.query(capital_country).areaId()
countries_and_relation_capitals.append([country[0], country[1], areaId, country[3]])
num_pharmacies = []
for country in countries_and_relation_capitals:
capitalId = country[2]
if capitalId is None:
continue
capitalId = capitalId
overpass = Overpass()
query = overpassQueryBuilder(area=capitalId, elementType=['node','way', 'relation'], selector='"amenity"="pharmacy"', out='count')
result = overpass.query(query, timeout=900)
num_pharmacies.append([result.countElements(), country])
def maximum_pharmacy(num_pharmacies):
maximum = num_pharmacies[0]
for country in num_pharmacies:
if maximum[0] < country[0]:
maximum = country
return maximum
max_pharmacy = maximum_pharmacy(num_pharmacies)
print(f'В {max_pharmacy[1][3]} столице {max_pharmacy[1][1]} содержится наибольшее количество аптек - {max_pharmacy[0]}')
bangladesh_id = max_pharmacy[1][0]
В Dhaka столице Bangladesh содержится наибольшее количество аптек - 1286
areaId = 3600000000 + bangladesh_id
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"amenity"="school"')
result = overpass.query(query, timeout=900)
points = []
for school in result.elements():
if school.lon() is None or school.lat() is None:
continue
points.append((school.lon(), school.lat()))
idx = rtree.index.Rtree()
for i, p in enumerate(points):
idx.insert(i, p + p, p)
p = (way[-1].x, way[-1].y)
five_nearest = list(idx.nearest(p, 5, objects='raw'))
orig_point = Point(way[-1])
lines = [LineString([orig_point, Point(dest)]) for dest in points if dest in five_nearest]
geod = Geod(ellps="WGS84")
lines_length = [geod.geometry_length(line) for line in lines]
print(f'{lines_length} расстояния от точки до 5 ближайших школ в Бангладеше')
[2284559.9538974287, 2266832.2552604834, 2281156.469616435, 2281167.7889800607, 2267737.838541723] расстояния от точки до 5 ближайших школ в Бангладеше
sum_length = 0
for i in lines_length:
sum_length += i
print(f'{sum_length} - сумма расстояний в метрах от точки до 5 ближайших школ в Бангладеше')
11381454.306296133 - сумма расстояний в метрах от точки до 5 ближайших школ в Бангладеше
import overpy
import shapely.geometry as geometry
from shapely.ops import linemerge, unary_union, polygonize
all_polygons = []
for area in areas:
query = """[out:json][timeout:25];rel(""" + str(area[0]) +""");out geom;>;out skel qt; """
api = overpy.Overpass()
result = api.query(query)
lss = [] #convert ways to linstrings
for ii_w, ways in enumerate(result.ways):
ls_coords = []
for node in ways.nodes:
ls_coords.append((node.lon, node.lat)) # create a list of node coordinates
lss.append(geometry.LineString(ls_coords)) # create a LineString from coords
merged = linemerge([*lss]) # merge LineStrings
borders = unary_union(merged) # linestrings to a MultiLineString
polygons = list(polygonize(borders))
all_polygons.append([area[0], area[1], geometry.MultiPolygon(polygons)])
south_and_west_for_countries = []
for multypolygon in all_polygons:
if multypolygon[1] == "France - Luxembourg" or multypolygon[1] == "France - Andorra":
continue
multy_poly = multypolygon[2]
x = multy_poly[0].exterior.coords.xy[0][0]
y = multy_poly[0].exterior.coords.xy[1][0]
south = [x, y]
north = [x, y]
for poly in multy_poly:
xs = poly.exterior.coords.xy[0]
ys = poly.exterior.coords.xy[1]
for x, y in zip(xs, ys):
if y < south[1]:
south = [x, y]
if y > north[1]:
north = [x, y]
south_and_west_for_countries.append([south, north, multypolygon[0], multypolygon[1]])
railway_station = []
for country in south_and_west_for_countries:
areaId = 3600000000 + country[2]
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node'], selector='"railway"="station"', out='geom')
result = overpass.query(query, timeout=900)
count = 0
if result.countElements() == 0:
continue
for node in result.nodes():
geod = Geod(ellps="WGS84")
south_line_length = geod.geometry_length(LineString([Point(country[0]), Point([node.lon(), node.lat()])]))
north_line_length = geod.geometry_length(LineString([Point(country[1]), Point([node.lon(), node.lat()])]))
if south_line_length < north_line_length:
count += 1
railway_station.append([count, country[2], country[3]])
bus_stop = []
for country in south_and_west_for_countries:
areaId = 3600000000 + country[2]
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node'], selector='"highway"="bus_stop"', out='geom')
result = overpass.query(query, timeout=900)
count = 0
if result.countElements() == 0:
continue
for node in result.nodes():
geod = Geod(ellps="WGS84")
south_line_length = geod.geometry_length(LineString([Point(country[0]), Point([node.lon(), node.lat()])]))
north_line_length = geod.geometry_length(LineString([Point(country[1]), Point([node.lon(), node.lat()])]))
if south_line_length < north_line_length:
count += 1
bus_stop.append([count, country[2], country[3]])
max_ratio = [railway_station[0][0] / bus_stop[0][0]]
for i in railway_station:
for j in bus_stop:
if i[2] == j[2]:
if j[0] == 0:
continue
ratio = i[0] / j[0]
if ratio > max_ratio[0]:
max_ratio = [ratio, i[0], j[0], i[1], i[2]]
print(f'В {max_ratio[4]} (id: {max_ratio[3]}) самое большое отношение числа железнодорожных станций к автобусным остановкам в южной части - {max_ratio[0]}')
cameroon_id = max_ratio[3]
В Cameroon (id: 192830) самое большое отношение числа железнодорожных станций к автобусным остановкам в южной части - 0.6181818181818182
areaId = 3600000000 + cameroon_id
overpass = Overpass()
query = overpassQueryBuilder(area=areaId, elementType=['node','way', 'relation'], selector='"amenity"="cinema"')
result = overpass.query(query, timeout=900)
points = []
for cinema in result.elements():
if cinema.lon() is None or cinema.lat() is None:
continue
points.append((cinema.lon(), cinema.lat()))
idx = rtree.index.Rtree()
for i, p in enumerate(points):
idx.insert(i, p + p, p)
p = (way[-1].x, way[-1].y)
five_nearest = list(idx.nearest(p, 5, objects='raw'))
orig_point = Point(way[-1])
lines = [LineString([orig_point, Point(dest)]) for dest in points if dest in five_nearest]
geod = Geod(ellps="WGS84")
lines_length = [geod.geometry_length(line) for line in lines]
sum_length = 0
for i in lines_length:
sum_length += i
print(f'{sum_length} - сумма расстояний в метрах от точки до 5 ближайших кинотеатров в Камеруне')
19612688.559784684 - сумма расстояний в метрах от точки до 5 ближайших кинотеатров в Камеруне
Добраться на автомобиле от входа в Central Park Нью-Йорка (со стороны 5th Avenue) до пересечения Water Street и Washington Street в Бруклине (откуда получаются лучшие фото Манхэттенского моста) довольно непросто - разумеется, из-за вечных пробок. Однако еще сложнее это сделать, проезжая мимо школ, где дети то и дело переходят дорогу в неположенном месте.
Вам необходимо построить описанный выше маршрут, избегая на своем пути школы. Визуализируйте данный маршрут (также добавив школы и недоступные для проезда участки дорог) при помощи Folium
Данные о расположении школ Нью-Йорка можно найти здесь
import json
import requests
import warnings
import folium
import pyproj
from shapely import geometry
from shapely.geometry import Point, LineString, Polygon, MultiPolygon
from openrouteservice import client
warnings.filterwarnings("ignore")
url = 'https://data.cityofnewyork.us/api/views/a3nt-yts4/rows.json?accessType=DOWNLOAD'
def CreateBufferPolygon(point_in, resolution=10, radius=10):
sr_wgs = pyproj.Proj(init='epsg:4326')
sr_utm = pyproj.Proj(init='epsg:32632') # WGS84 UTM32N
point_in_proj = pyproj.transform(sr_wgs, sr_utm, *point_in) # unpack list to arguments
point_buffer_proj = Point(point_in_proj).buffer(radius, resolution=resolution) # 10 m buffer
# Iterate over all points in buffer and build polygon
poly_wgs = []
for point in point_buffer_proj.exterior.coords:
poly_wgs.append(pyproj.transform(sr_utm, sr_wgs, *point)) # Transform back to WGS84
return poly_wgs
# Set up the fundamentals
api_key = '5b3ce3597851110001cf624857c34d71aefe42ed95ec0fb54e33a8ea' # Individual api key
clnt = client.Client(key=api_key) # Create client with api key
school_json = requests.get(url).json() # Get data as JSON
map_params = {'location':([40.703165, -73.987598]),
'zoom_start': 5}
map1 = folium.Map(**map_params)
# Populate a school buffer polygon list
schools_poly = []
for school_data in school_json['data']:
school_coords = school_data[8].split()
school_coords = [float(school_coords[1][1:]), float(school_coords[2][:-1])]
folium.features.Marker(list(reversed(school_coords)),
popup='School point<br>{0}'.format(school_coords)).add_to(map1)
# Create buffer polygons around construction sites with 150 m radius and low resolution
school_poly_coords = CreateBufferPolygon(school_coords,
resolution=2, # low resolution to keep polygons lean
radius=150)
schools_poly.append(school_poly_coords)
school_poly_coords = [(y,x) for x,y in school_poly_coords] # Reverse coords for folium/Leaflet
folium.vector_layers.Polygon(locations=school_poly_coords,
color='#ffd699',
fill_color='#ffd699',
fill_opacity=0.2,
weight=3).add_to(map1)
# GeoJSON style function
def style_function(color):
return lambda feature: dict(color=color,
weight=3,
opacity=0.5)
# Create new map to start from scratch
map_params = {'location':([40.744696, -73.953993]),
'zoom_start': 12}
map2 = folium.Map(**map_params)
request_params = {'coordinates': [[-73.963793, 40.777111],
[-73.989596, 40.703213]],
'format_out': 'geojson',
'profile': 'driving-car',
'preference': 'shortest',
'instructions': 'false',}
route_normal = clnt.directions(**request_params)
folium.features.GeoJson(data=route_normal,
name='Route with school',
style_function=style_function('#FF0000'),
overlay=True).add_to(map2)
school_buffer_poly = []
for school_poly in schools_poly:
poly = Polygon(school_poly)
school_buffer_poly.append(poly)
# Add the school polygons to the request parameters
request_params['options'] = {'avoid_polygons': geometry.mapping(MultiPolygon(school_buffer_poly))}
route_detour = clnt.directions(**request_params)
#Buffer route with 0.009 degrees (really, just too lazy to project again...)
route_buffer = LineString(route_detour['features'][0]['geometry']['coordinates']).buffer(0.009)
folium.features.GeoJson(data=geometry.mapping(route_buffer),
name='Route Buffer',
style_function=style_function('#FFFF00'),
overlay=True).add_to(map2)
for school_poly in schools_poly:
poly = Polygon(school_poly)
if route_buffer.intersects(poly):
school_poly_coords = [(y,x) for x,y in school_poly] # Reverse coords for folium/Leaflet
folium.vector_layers.Polygon(locations=school_poly_coords,
color='#ffd699',
fill_color='#ffd699',
fill_opacity=0.2,
weight=3).add_to(map2)
folium.features.Marker(list(reversed(poly.centroid.coords[0]))).add_to(map2)
folium.features.GeoJson(data=route_detour,
name='Route without schools',
style_function=style_function('#00FF00'),
overlay=True).add_to(map2)
<folium.features.GeoJson at 0x7fe054a59dd0>
Синими маркерами отмечены школы в Нью-Йорке, которые потенциально могли бы встретиться; оранжевыми многоугольниками отмечены недоступные для проезда участки дорог. Зеленой линией отмечен маршрут, на пути которого нет школ; красной линией отмечен маршрут, в котором НЕ учитывается расположение школ
Все школы не были указаны на карте по причине их большого количества. Jupyter Notebook начинает подлагивать, если честно указывать все школы в Нью-Йорке
map2
В следующих строках происходит нанесение всех школ на карту. Получается в точности такая же карта, что и map2, за исключением того, что на карту нанесены все школы в Нью-Йорке. Получившуюся карту (map3) лучше не выводить, так как все начинает сильно лагать
# GeoJSON style function
def style_function(color):
return lambda feature: dict(color=color,
weight=3,
opacity=0.5)
# Create new map to start from scratch
map_params = {'location':([40.744696, -73.953993]),
'zoom_start': 12}
map3 = folium.Map(**map_params)
request_params = {'coordinates': [[-73.963793, 40.777111],
[-73.989596, 40.703213]],
'format_out': 'geojson',
'profile': 'driving-car',
'preference': 'shortest',
'instructions': 'false',}
route_normal = clnt.directions(**request_params)
folium.features.GeoJson(data=route_normal,
name='Route without construction sites',
style_function=style_function('#FF0000'),
overlay=True).add_to(map3)
# Plot which schools fall into the buffer Polygon
school_buffer_poly = []
for school_poly in schools_poly:
poly = Polygon(school_poly)
folium.features.Marker(list(reversed(poly.centroid.coords[0]))).add_to(map3)
school_poly_coords = [(y,x) for x,y in school_poly] # Reverse coords for folium/Leaflet
folium.vector_layers.Polygon(locations=school_poly_coords,
color='#ffd699',
fill_color='#ffd699',
fill_opacity=0.2,
weight=3).add_to(map3)
school_buffer_poly.append(poly)
# Add the school polygons to the request parameters
request_params['options'] = {'avoid_polygons': geometry.mapping(MultiPolygon(school_buffer_poly))}
route_detour = clnt.directions(**request_params)
folium.features.GeoJson(data=route_detour,
name='Route with construction sites',
style_function=style_function('#00FF00'),
overlay=True).add_to(map3)
# map3.add_child(folium.map.LayerControl())
<folium.features.GeoJson at 0x7fe04c482fd0>
Синими маркерами отмечены школы в Нью-Йорке; оранжевыми многоугольниками отмечены недоступные для проезда участки дорог. Зеленой линией отмечен маршрут, на пути которого нет школ; красной линией отмечен маршрут, в котором НЕ учитывается расположение школ
route_buffer = LineString(route_detour['features'][0]['geometry']['coordinates']).buffer(0.009)
folium.features.GeoJson(data=geometry.mapping(route_buffer),
name='Route Buffer',
style_function=style_function('#FFFF00'),
overlay=True).add_to(map3)
<folium.features.GeoJson at 0x7fe0326cb090>
map3